#include <iostream>


#include "Matrix.h"


class SLE 	    //system of linear equations n X n
{
    SquareMatrix A;
    Matrix B;

	static Matrix jacobi(const SquareMatrix& A,const Matrix& B);
	static Matrix gaussSiedel(const SquareMatrix& A,const Matrix& B);
	
	bool securePivot(long row);
	long findGoodRow(long row) const;
	bool makeOne(long row);
	bool makeZero(long row);
	
 public:
    
    SLE(){}
    SLE(const char* matA,const char* matB);
    SLE(SquareMatrix A,Matrix B); 
    SLE(const SLE& another);
    ~SLE();
    SLE& operator=(const SLE& another);
    
    void useIterativeMethods();
    Matrix gaussianElimination();
	
	
};




// ########  SLE method defifnitions  #######################

SLE::SLE(const char * matA,const char* matB )
{
	A= SquareMatrix(matA);
	B= Matrix(matB);

}
	    
	    
SLE::SLE(SquareMatrix A,Matrix B )
{
	this->A=A;
	this->B=B;

}

SLE::SLE(const SLE& another)
{
	A=another.A;
	B=another.B;

}


SLE::~SLE()
{

}

SLE& SLE::operator=(const SLE& another)
{
	A=another.A;
	B=another.B;

	return (*this);
}



void SLE::useIterativeMethods()
{
	if(!A.makeDiagonallyDominant(B))
	{
	    cout <<"\nJacobi and Gauss-Siedel methods cannot be applied as matrix A cannot be made diagonally dominant.\n";
		return;
	}
	
	jacobi(A,B).disp();
	gaussSiedel(A,B).disp();

}


Matrix SLE::jacobi(const SquareMatrix& A,const Matrix& B)
{
		
    SquareMatrix T=A;
	Matrix c=B;
	
	double p;
	for(long i=0;i<T.numRows();++i)
	{
	    p=T[i][i];
	    T.multRow(i,-1/p);
		c.multRow(i,1/p);
		T[i][i]=0;
	}
	
	Matrix x=c;
	Matrix x1=x;
	Matrix solution(c);
	
	long iterations=0;
	do
	{
	    iterations++;
	
	    x=x1;
	    x1=T*x +c;
		
		solution=solution.concat(x1);
		
	}while(!x1.nearlyEqualTo(x));
	
	cout <<"\n number of iterations for Jacobi : " <<iterations <<"\n";
	return solution;
}



Matrix SLE::gaussSiedel(const SquareMatrix& A,const Matrix& B)
{
    SquareMatrix T=A;
	Matrix c=B;
	
	double p;
	for(long i=0;i<T.numRows();++i)
	{
	    p=T[i][i];
	    T.multRow(i,-1/p);
		c.multRow(i,1/p);
		T[i][i]=0;
	}
	
	Matrix x=c;
	Matrix x1(x);
	Matrix solution(c);
	long iterations=0;
		
	do
	{
	    x=x1;		
		iterations++;
		
	    for(long i=0;i<T.numRows();++i)
		{
            x1[i][0]=(Matrix(T[i]) * x1)[0][0] + c[i][0];
		}
		
		solution=solution.concat(x1);
		
	}while(!x1.nearlyEqualTo(x));
			
	cout <<"\n number of iterations for Gauss-Siedel : " <<iterations <<"\n";
	
	return solution;
    
	
	/*  the code when i did not have Matrix(Vector) constructor
	do
	{
	    x=x1;		
		iterations++;
		
	    // processing without Matrix(Vector) constructor
	    for(long j=0;j<T.numRows();++j)
		{
		    tempx=0;
		    for(long k=0;k<T.numCols();++k)
		    {
                tempx=tempx + T[j][k] * x1[k][0] ;
			}
			
			x1[j][0]=tempx+c[j][0];
		}
		
		
	}while(!x1.nearlyEqualTo(x));
	*/
	
	
}


Matrix SLE::gaussianElimination()
{
	
	long till=A.numRows();
	for(long i=0;i<till;++i)
	{
		if(Vector::nearlyEqual(A[i][i],0.0) and !securePivot(i))
		{
			cout <<"cannot apply Gaussian Elimination as the pivot cannot be made non-zero\n";
			return false;
		}
	
		makeOne(i);
		makeZero(i);
	}

	cout.precision(10);
	for(long i=0;i<till;++i)
	{
		cout <<B[i][0] <<"  ";
	}
	
	cout <<"\n";
	
	return B;
}


bool SLE::securePivot(long row)
{
	long goodRow;
	if((goodRow=findGoodRow(row))!=-1)
	{
		A.swapRows(row,goodRow);
		B.swapRows(row,goodRow);
	}
	else
	{
		return false;
	}

	return true;
}


long SLE::findGoodRow(long row) const
{
	long till=A.numRows();
	
	for(long i=row+1;i<till;++i)
	{
		if(!Vector::nearlyEqual(A[i][row],0.0))
		{
			return i;
		}
		
	}
	
	return -1;	

}


bool SLE::makeOne(long row)
{
	double multiplier=1/A[row][row];
	
	if(multiplier==1)
		return true;
		
	A.multRow(row,multiplier);
	B.multRow(row,multiplier);
	
	return true;

}


bool SLE::makeZero(long row)
{
	long till=A.numRows();
	
	for(long i=0; i<till;++i)
	{
		if(i==row)
			continue;
			
		double multiplier=A[i][row]*(-1);	
		A.addMult(i,multiplier,row);
		B.addMult(i,multiplier,row);
	}
		
	return true;
		
}



//  ###########   end of method definitions  ###################


/*
int main()
{
    SLE sle("A5","B5");
	
	sle.gaussianElimination();
	//sle.useIterativeMethods();	
	
	return 0;
}
*/

